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Abstract 

Excited-state calculations, notably for quasiparticle band structures, are nowadays 
routinely performed within the GW approximation for the electronic self-energy. 
Nevertheless, certain numerical approximations and simplifications are still em- 
ployed in practice to make the computations feasible. An important aspect for 
periodic systems is the proper treatment of the singularity of the screened Coulomb 
interaction in reciprocal space, which results from the slow 1/r decay in real space. 
This must be done without introducing artificial interactions between the quasi- 
particles and their periodic images in repeated cells, which occur when integrals of 
the screened Coulomb interaction are discretised in reciprocal space. An adequate 
treatment of both aspects is crucial for a numerically stable computation of the 
self-energy. In this article we build on existing schemes for isotropic screening and 
present an extension for anisotropic systems. We also show how the contributions 
to the dielectric function arising from the non-local part of the pseudopotentials can 
be computed efficiently. These improvements are crucial for obtaining a fast conver- 
gence with respect to the number of points used for the Brillouin zone integration 
and prove to be essential to make GW calculations for strongly anisotropic systems, 
such as slabs or multilayers, efficient. 
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1 Introduction 



For describing quasiparticle excitations, as measured in direct and inverse pho- 
toemission, many-body perturbation theory in the GW approximation [1] has 
developed into the method of choice for weakly correlated solids and their sur- 
faces. In particular, the GW approximation describes the quasiparticle band 
structures and band gaps for a large variety of semiconductors in good agree- 
ment with experimental results [2,3]. For more correlated systems, it becomes 
necessary to go beyond GW, but these schemes often include the GW self- 
energy diagrams as lowest order [4,5,6]. Similarly, the Bethe-Salpeter approach 
to electron-hole excitations, as probed in optical absorption or electron energy- 
loss spectroscopy, builds on the GW self-energy [7]. 

In the GW approximation the frequency-dependent, non-local self-energy E 
that connects the independent-particle Green function Go with the interacting 
one G is given by S = iGW, where W is the dynamically screened Coulomb in- 
teraction. The independent-particle starting point is typically chosen to be the 
Green function of a Kohn-Sham density-functional theory (DFT) calculation. 
Despite significant methodological progress [8,9,10,11,12,13,14] GW calcula- 
tions are still computationally demanding, and their application is limited 
to relatively small system sizes. So far all implementations employ a number 
of additional simplifications to reduce the computational cost. Some of these 
are motivated by physical considerations, such as plasmon-pole models [15] or 
model dielectric functions [16,17,18,19], while others appear as purely math- 
ematical "tricks" to improve the numerical stability or efficiency Often, the 
validity and usefulness of a specific approach depends on the physical sys- 
tem under consideration. An important aspect in every GW implementation 
that uses reciprocal space is the treatment of the singularity at k — > in 
the bare and screened interaction for non-metallic systems. This singularity 
is integrable, and many different schemes have been developed for these k- 
space integrals [11,15,20,21,22,23]. The general idea is to describe the singular 
part by a model function that can be handled analytically, so that the re- 
mainder is sufficiently smooth for a numerical treatment. In physical terms, 
the k — > behaviour determines the long-range part of the interaction. The 
discretisation of the reciprocal-space integrals can then be interpreted as in- 
troducing an artificial supercell periodicity in real space. The periodic images 
of the quasiparticles give rise to an infinite self-interaction due to the l/|r| 
tail of the interaction. Any scheme that integrates the singularity analytically 
is therefore strictly equivalent to modifying the long-range behaviour of the 
interaction such that a quasiparticle does not interact with its periodic im- 
ages or that the interaction decays faster than l/|r| 2 . In other words, different 
integration schemes correspond to particular modifications of the long-range 
tail and vice versa. 
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In most schemes screening is assumed to be isotropic at the length scale given 
by the inverse of the smallest non-zero k-vector, which is appropiate for most 
bulk materials but may fail in systems with an appreciable anisotropy, such as 
superlattices or layered materials, as well as in supercell approaches for low- 
dimensional materials like clusters, molecules, nanowires, films, or surfaces. 
An obvious way to avoid these spurious interactions for systems with broken 
translational symmetry is to abandon the concept of periodic boundary condi- 
tions altogether in the relevant directions and perform the calculation entirely 
in real space. For semi-infinite jellium surfaces such a GW embedding scheme 
has been successfully implemented [24,25]. Its extension to realistic surfaces, 
however, is computationally still too expensive. A real-space implementation 
for finite systems has also been reported [13,14], but its applicability to sys- 
tems with periodicity in one or more directions remains to be shown. 



Staying with the repeated-cell approach, we will show in the following how 
it is possible to incorporate the anisotropy in the treatment of the singu- 
larity in the GW space-time method [20]. In addition to the equations that 
we have implemented we derive exact expressions that allow us to discuss 
other algorithms in comparison. Furthermore, we will show that the proposed 
modifications considerably improve the convergence behaviour with respect 
to the number of k-points, which is the natural parameter associated with 
the singularity treatment. In practice the GW approximation is often applied 
non-self-consistently by constructing the screened interaction as well as the 
self-energy from the independent-particle Green function Gq. However, the 
behaviour of the anisotropy discussed in this article applies equally to the 
fully self-consistent GW approach. For simplicty we will therefore focus on 
the non-self-consistent case and indicate differences whenever they apply. In 
the interest of readability we will also refrain from introducing different sym- 
bols to distinguish between the self and non-self-consistent case. 



This paper is organised as follows. In Section 2 we briefly describe the com- 
putation of the self-energy in the space-time method. In Section 3 we explain 
how the anisotropy is accounted for; the detailed derivation of the well-known 
anisotropic equations is presented in Appendix B for completeness. In Sec- 
tion 4 we demonstrate the improved k-point convergence behaviour result- 
ing from our modifications before summarising our results in Section 5. In 
Appendix A we have collected the spherical-harmonics expansion of several 
vector quantities that appear in our derivations. Finally, an efficient imple- 
mentation of the contribution from Kleinman-Bylander-type non-local pseu- 
dopotentials [26], which enter the expressions for the anisotropy, is presented 
in Appendix C. Unless otherwise indicated, we use Hartree atomic units. 
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2 Outline of the GW space-time method 



The GW space-time method has been presented in detail elsewhere [20,27]. 
We will therefore only summarise the steps to construct the self-energy from 
the output of a preceding DFT calculation. Assuming a non-magnetic systems 
for simplicity (the extension to a spin-dependent Green function is straight- 
forward), the computational steps are: 

(1) Construction of the non- interacting Green function G in real space and 
imaginary time from the Kohn-Sham eigenfunctions <^ n k and eigenvalues 
e„k (the Fermi level is set as the energy zero) 

!occ 
EvWrR k (r')e— , r < 0, 
(1) 
- E </WrR k (r')e-^ , r > 0, 
n 

where Q denotes the unit-cell volume and the integral over k runs over 
the first Brillouin zone, 

(2) formation of the irreducible polarisability P in the random-phase approx- 
imation in real space and imaginary time 

P(r, r'; ir) = -2iG(r, r'; ir)G(r', r; -ir) , (2) 

(3) Fourier transformation of P to reciprocal space 

P GG ,(k, ir) = ^J d 3 r J dV P(r, r'; ir ) e - i (k+G). r+ i(k + G').r' (3) 

and to imaginary frequency, 

(4) construction of the symmetrised dielectric matrix in reciprocal space 

An 

e gg , (k, iu) = 5 GG > - | k+G || k+G /| P GG' (k, iu) , (4) 

(5) inversion of the symmetrised dielectric matrix for each k-point and each 
imaginary frequency, 

(6) calculation of the screened Coulomb interaction in reciprocal space 

Att 

W GG , (k, iu) = | k + G || k+G /| g GG' ( k > i^) , (5) 

(7) Fourier transformation of W to imaginary time and to real space 

W(r, r'; ir) = J— [d 3 k ^ W GG ,(k, ir)e Kk + G). r -i( k+ G0. r ' > (6) 

y Z7r ) B J Z g,c 
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(8) computation of the self-energy in real space and imaginary time 



S(r, r'; ir) = iG(r, r'; W)W(r, r'; ir) . 



(7) 



The Coulomb singularity appears explicitly for G = or G' = as k — > 
in steps 4 and 6; in the actual implementation, however, it is treated in 
steps 5 and 7 for numerical reasons. The anisotropy enters the scheme naturally 
through the construction of the dielectric matrix and must be taken fully into 
account in the screened interaction. 

The quasiparticle energies are obtained by computing the matrix elements of 
the self-energy (v?nk|£(ir)|v? n k) on the imaginary time axis, which are then 
Fourier-transformed to imaginary frequency and analytically continued to the 
real frequency axis. Approximating the quasiparticle by the DFT Kohn-Sham 
wavefunctions finally gives the quasiparticle energies e^ k as solutions of the 
quasiparticle equation 



where V xc is the exchange-correlation potential used in the underlying DFT 
calculation. The details of the analytic continuation and the solution of Equa- 
tion (8) have been described elsewhere [20]. Self-consistency in GW would be 
achieved by entering step 2 with a new Green function obtained from solving 
Dyson's equation G = G + CqSG after step 8 and iterating steps 2-8. 



3 Treatment of the anisotropy 

3.1 Anisotropy in the screened interaction 

As shown in Appendix B the head (G = G' = 0) and the wings (G = or 
G' = 0) of the dielectric matrix close to the T-point, i.e., for k — > 0, depend 
on the direction in which this limit is taken. We denote this dependence by the 
spatial angle f2k, and the corresponding normalised direction vector by k. For 
simplicity, the imaginary frequency argument iu is omitted in the following. 
The directional dependence at the T-point is present in the whole inverse 
dielectric matrix, i.e., head, wings, and body. By block- wise inversion [28] it is 
easily shown that the head of the inverse symmetrised dielectric matrix takes 
the form (cf. Appendix B) 



(8) 



£ Oo(^k) 



1 



(9) 



k T Lk 
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where the matrix L is the macroscopic dielectric tensor. We note that in most 
other implementations, in which this anisotropy has been considered, such as 
[22,23], but not [21], the right-hand side of Equation (9) has been replaced by 
the expression k T L _1 k without formal justification. 



Correspondingly, the wings can be expressed as 
^(nk) = -e5b 1 ("k)fk.S(G)" . 



(10) 



The vector S(G) is defined in Equation (B.15). For the remainder of this 
article we will restrict the discussion of the wings to the expression for G' = 
since the case G = is trivially obtained from the symmetry relation 



^OG'(^k) - kc'o(^k)]* • 



Finally, the body is given by 

^GG'(^k) = ^GG' + ^00 (^k) 



k-SfG' 



k-SfG"! 



:i2) 



where B denotes the body of the symmetrised dielectric matrix as defined in 
Equation (B.ll). Expressions (9)-(12) imply that each element of the inverse 
dielectric matrix is, in general, not analytic at k = 0, i.e., it does not have a 
unique limit for k — > 0. This was already recognised more than 30 years ago 
by Pick et al. [28], but the treatment of this non-analytic behaviour in GW 
implementations has not been discussed widely in the literature. Hybertsen 
and Louie address the problem in the appendix of [15], but in the actual 
calculation they neglect the non-analytic part, arguing that the error can be 
made negligibly small with a sufficiently high number of k-points. So far the 
anisotropy has only been considered for the head element in connection with 
the treatment of the Coulomb singularity [21,22,23]. We will return to this 
point later in Section 3.5. 

Combining Equations (9) to (12) with Equation (5), we obtain the screened 
interaction for k — > 



W 00 (k) 
W G0 (k) 
W GG ,(k) 



47T 

|kj2 e °° 

47T 



|k||G| 

47T 



(fik), 



e5}(n k ) k-S(G) 



IGIIG' 



7T (Sgg' + £ ~°o(^k) [ k • S ( G )] [ k ' S(G')] *) 



(13) 
(14) 
(15) 



6 



The presence of the singularity at G = G' = (Equation (13)) necessi- 
tates a special numerical treatment. In the space-time method, this problem 
is solved by splitting off a long-range part W lr with the appropiate k — > 
behaviour, which is chosen such that its Fourier transform can be computed 
semi-analytically as described in detail in the following section. The remaining 
short-range part W sr = W — W h can then safely be treated numerically since 
it is no longer singular. 

The next step in the space-time method is the Fourier transformation of W to 
real space (Equation (6)); in reciprocal-space algorithms it is the construction 
of the matrix elements of E. Both approaches involve an integration over 
the Brillouin zone after multiplication by an analytic function a GG /(k), and 
it is in this integration that the anisotropy must be taken into account. In 
practice, these integrals are usually discretised, which we express formally by 
partitioning the Brillouin zone into subzones Zj with volume Vf 

Jd 3 k W GG ,(k)a GG ,(k) = £ /d 3 fc W GG ,(k)a GG ,(k) • (16) 

BZ * Zj 

We denote the subzone that contains the T-point by Z-p and assume that it 
has inversion symmetry about k = 0. While the subzone integrals for i ^ T 
can be approximated by 

J d 3 k W GG > (k)a GG / (k) » WGG'(k*)GGG'(kO , (17) 

z, 

where k, is a representative point for the subzone Zj - usually its centre - the 
integrals over Zr require a special treatment due to the non-analyticity of W 

at r. 

In principle, even the singularity for G = G' = can be treated in this way. 
Since it is illuminating to discuss existing isotropic and anisotropic singularity 
integration schemes in terms of approximations to an exact expression, we 
present the corresponding equations for reciprocal space algorithms in Section 
3.5. In the space-time method, on the other hand, the separation into W lr and 
W ST is more efficient than the direct approach, because the analytic function 
°GG'(k) for the Fourier transformation 

a GG ,(k) = e i ( k+G )- r - i ( k+G '> r ' (18) 

depends on r and r' and would thus require the computation of the integrals 
over Z r in Equation (16) for every r and r'. 
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Since the long-range part of the interaction yields a significant contribution 
to the quasiparticle energies, an accurate treatment of its anisotropy is very 
important and is therefore described first in Section 3.2. In Section 3.3 we 
show that the apparent l/|k| singularity of the wings does not cause numerical 
problems for the integrals over Zr, and in Section 3.4 the computation of the 
integrals for the body is presented. 



3.2 Treatment of the head 



In the space-time method the head of the inverse dielectric matrix is used to 
define the long-range part of the screened interaction. For this purpose, we 
extend Equation (13) to G = G' ^ and define the long-range part for all k 
in the Brillouin zone as 

H "-' (k » = ( k + TO + G ) fe °- • (19) 



For numerical reasons we subtract the long-range part at the level of the in- 
verse dielectric matrix after applying the body corrections described in Section 
3.4 for k = 0, and compute W sr from this modified entity according to 



— l,sr/i\ _ — 1 /. n |k + G| 

£ GG ,W:-£ GG , V „, (k+G) T L(k + G) ' 



£(k) :=s~ GG ,(k) - —±^-l——S GGf , (20) 



= • < 21 > 

By expanding the angular dependence of W lr into spherical harmonics (cf. 
Appendix A) 

00 ' Aix 

WgoOO = E E H lm 5 GG/ Y lm (n k+G ) , (22) 

Z=0 m=—l I I 



the Fourier transformation of W can be performed analytically. Only even / 
contribute to the sum because the coefficients H\ m vanish for odd I. Making use 
of the expansion of a plane wave [29] in spherical harmonics Y% m and spherical 
Bessel functions j t , 

oo I 

e ik ' r = 4ttE E i l Ji(kr)Y lm (^)Y; m (n k ) (23) 

1=0 rn=-l 
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we arrive at 



OO t 1 

^ lr (r,r') = E E cjH^YUn^,) — 



;=0 m=-/ 



(24) 



The coefficients q for even / are defined as 



oo 



Q = - 




(25) 



71 J 




with n\\ = n(n — 2)(n— 4) • • •. In practice we truncate the sum in Equation (24) 
at finite / = Z max (cf. Section 4.1). 

For numerical convenience £ is split into a static exchange part S x = iGv 
and a frequency-dependent correlation part S c = iG(W — v) in the space- 
time method [20]. This is achieved by subtracting the unscreened Coulomb 
interaction v from W lr in its angular expansion (22), i.e., is subtracted 

from Hqq for each imaginary frequency. Furthermore, the transformation from 
imaginary frequency to imaginary time is performed on the expansion coeffi- 
cients Hi m {iuj) directly, and we obtain (W h — v) according to Equation (24) 
with the expansion coefficients in imaginary time Hi m {ir). 

A proper treatment of the anisotropy in the long-range part of the screened 
interaction is crucial to obtain converged results. This is easily illustrated in 
the space-time method: For non-local operators like W or E the density of 
the k-point sampling determines the range of the non-locality in real space. 
For example, a 4x4x4 k-grid corresponds to a maximum non-locality range 
or interaction cell of 4 real-space unit cells in each dimension. If parts of the 
long-range interaction remain in W sr for small but finite k, the tails of W ST 
extend over the boundary of the interaction cell and will be folded back in 
the numerical Fourier transformation when applying the periodic boundary 
conditions. Since the size of the interaction cell is determined by the k-point 
sampling, an inadequate treatment of the long-range part would result in an 
unsatisfactory k-convergence behaviour. 

3. 3 Treatment of the wings 

The wings are antisymmetric with respect to k, i.e., 



W^Go(k) 



W^Go(-k) . 



(26) 



Hence we can write the T-point contribution as 
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|d 3 A;^ G o(k)a G o(k) 

\(f d 3 kW G0 (k)a G0 (k)+ [ d 3 kW G0 (-k)a G0 (-k) 



47T 

IgT 



d^eoo 1 ^) k-S(G) k-V k a G0 (fi k )| k=0 + O(|k| 2 ) 



(27) 



Zr 



where we have made use of the Taylor expansion of the analytic function 
dGo(k). The important benefit of this reformulation is that no term in the 
integrand is singular, and therefore we do not expect any numerical difficulties 
close to the T-point. This applies equally to the Fourier transformation in 
the space-time method as well as to the construction of the self-energy in 
reciprocal-space approaches. 

In practice we neglect the wing contributions from the T-point in the Fourier 
transformation, because evaluating Equation (27) for each r and r' would be 
computationally very demanding. The associated error scales as Vr- It is thus 
automatically controlled by the standard k-point convergence tests, since Vr 
is inversely proportional to the number of points in our regular k-point grid. 
Test calculations indicate that the overall convergence behaviour shows no 
significant improvement for a more sophisticated treatment of the wings. We 
note that in contrast to head and body the isotropic average for the wing 
contribution vanishes, i.e., no analytic contribution has to be considered when 
the non-analytic part is neglected. 



3.4 Treatment of the body 



For |k| <C |G| we see from Equation (18) that we can approximate a GG '(k) m 
a GG /(0). Similar considerations apply to the analytic function in reciprocal- 
space approaches. We can then express 

/ d 3 A; WW(k)a GG ,(k) « a GG ,(0) / d 3 A; ^^yp ( 2 §) 



as a discretised contribution analogous to Equation (17) with an averaged 
anisotropy 



£ GG' 



0) ;= -L |d 3 A;e GG ,(a). (29) 

r z r 
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We will now show how this average can be computed efficiently. To this end 
we rewrite the integral in polar coordinates as 

"■max (n k ) 

£gg'(0) = ydfi k w(^ k )£ G G'(^k) with w(tt k ) = — J k 2 dk , (30) 



where k max (Qk) is the distance from the centre to the surface of Zr in the 
direction of f2k, and w(fik) ac ts as an angular weight function that takes the 
shape of the T-zone element into account and may be subject to additional 
approximations, e.g., for spherical averages it is simply a constant u>(£V) = 

1/47T. 

Inserting Equation (12), we obtain a very simple expression 

oo I 1 1 

4'(0)=Bgg- + EE E E H lm s m ,(G)s;AG') 

1=0 m=—l m'=— 1 m"=— 1 

x jdQ k Y lm (n k )Y lml (Q k )Y* m „(Q k ) (31) 

when expanding [w(Q k )£oo(^k)] as well as k • S(G) in spherical harmonics as 
described in Appendix A. The angular integrals in Equation (31) are nothing 
but the Clebsch-Gordan coefficients for the spherical harmonics (Im lm'jlm"). 
From the properties of the Clebsch-Gordan coefficients [30] it follows that only 
a small number of non-zero terms contribute to Equation (31), namely I = 
with 3 terms and I = 2 with 9 terms. We refer to these 12 terms as "body 
corrections". In the original space-time implementation [20], Equation (12) 
was averaged over the Cartesian directions, which is equivalent to including 
only the / = terms with an approximate coefficient H 00 , calculated with the 
spherical weight function w(Q) = 1/Att and a 3-point integration. We note 
that the exact procedure includes only 9 more terms with 1 = 2. 



3. 5 Treatment of the anisotropy in reciprocal- space approaches 



For completeness, we present here a simple recipe to take the anisotropy into 
account in reciprocal-space approaches. The analytic function that appears 
in the computation of, e.g., the self-energy matrix element (0 nq |E(o/)|(/> nq ) is 
given by an expression of the form 

a GG ,(k) = [ d 8 gX;W(q.k)W(q.k)%w'.Wk) > (32) 

J m 



11 



where Mg n (q,k) = (0 nq _ k |e i(k+G)r |0 nq ) and F(u, u', e m(q _ k )) contains pref- 
actors and frequency integrals [11]. For body and wings the approach outlined 
in the previous sections also applies for this analytic function. For the head 
element, on the other hand, the integral to compute is 



d>* ^§ , (33, 



Z r 



where aoo(k) is analytic at k = 0. Assuming that aoo(k) is non-zero and varies 
sufficiently slowly with k, we set aoo(k) ~ aoo(0) and evaluate the remaining 
integral in polar coordinates 

. . "-max 

(«k) 

Z r 



In accordance with the computation of H 00 as described in Appendix A, the 
angular integral can be computed numerically on an appropriate angular grid 
with the angular weight function 

(«k) 

K(n k ) = J k 2 dk — = A; max (fi k ) . (35) 
o 



Alternatively, but formally equivalent, K(fl^) can be expanded in spherical 
harmonics with coefficients Ki m . The integral then becomes 



oo I 



jd 3 k ^ = aoo(0) E ± K; m H lm . (36) 

Z r '=0 m=-l 



Equation (36) is conducive for a discussion of different anisotropy and the 
Coulomb singularity treatments. In all isotropic approximations the sum over 
/ and m is restricted to the I — 0, m — term. In the "spherical" approximation 
[31] used in the early days of modern GW calculations [15,32] K 00 is further 
replaced by a shape-independent term 

1 / r \i/3 
^ h = ^r- {J dn ^ax(^)) • (37) 



The numerical computation of H 00 is restricted to 3 angular points if the 
average over the Cartesian directions is taken, which might introduce an addi- 
tional inaccuracy for anisotropic systems. In the improved integration scheme 
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used by Pulci et al. [22] as well as in the integration scheme by Wenzien et al. 
[23] , the head of the inverse dielectric matrix is written as a tensor and hence 
includes also the / = 2 terms (cf. Appendix A). The tensor itself is chosen to 
reproduce the correct value in the main directions, so Hi m is effectively de- 
termined from three independent points only. In the scheme that we propose 
here only the choice of the angular grid determines the accuracy of the sum 
in Equation (36). It should be comparable to the scheme proposed by Hott 
[21], which formally includes all terms and also involves a numerical integra- 
tion, the details of which are unfortunately not specified in [21]. Similarly, 
the offset-r-point method described in [11] in principle allows to capture the 
anisotropy to arbitrary precision. However, the accuracy with which this is 
achieved in practice depends on the choice of k-points. 



4 Results 

The equations presented in the previous sections were implemented into the 
gwst code [20,27]. The test system was chosen to be a periodically repeated 
4- layer Si (001) slab saturated with hydrogen and a vacuum separation be- 
tween the Si surface atoms equivalent to 4 layers. 15 points per half-axis were 
used both for the imaginary time and frequency Gauss-Legendre grids at a 
maximum numerical range of 6 atomic units. Convergence in the plane-wave 
cutoff is achieved for 7 Hartree, and unoccupied states up to 5 Hartree above 
the Fermi energy were included (610 bands). Head and wings were computed 
separately with a k-point grid of 14 x 14 x 1 and 120 bands. These parameters 
are sufficient to obtain quasiparticle energies converged to within 0.05 eV. 

Taking the body corrections from Equation (31) into account changes the 
quasiparticle energies of our test system only little compared to neglecting 
them completely and setting £gg'(0) = B^ G ,. The magnitude of the correc- 
tions depends on the weight of the T-point and hence on the k-point sampling, 
amounting to ~10meV for a 3 x 3 x 1 sampling and 1-2 meV for 8 x 8 x 1. 
While these corrections are small compared to the accuracy of our test calcu- 
lation, they might be larger for systems with stronger local-field effects. Since 
the additional computational effort is small we always include them. 

The repeated-slab arrangement considered here is a hypothetical system, that, 
at least until now, cannot be prepared in experiment. However, since it is fun- 
damentally a strongly anisotropic system, it provides an ideal test case for our 
modifications. It also provides the possibility to tune the degree of anisotropy 
by varying the slab thickness and the separation. If the slab separation were 
increased to infinity, the limit of an isolated slab would be recovered. If addi- 
tionally the slab thickness were taken to infinity the limit of a hydrogenated 
silicon surface would be reached. For the present choice of slab thickness and 
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Table 1 

Dependence of the quasiparticle energies (in eV relative to the valence-band maxi- 
mum) on the maximum angular momentum I in Equation (24). The k-point sam- 
pling for the data presented here is 4 x 4 x 4. Other samplings show the same 
behaviour. 



^max 





2 


4 


6 


lowest valence state 


-10.448 


-10.338 


-10.336 


-10.336 


lowest conduction state 


3.560 


3.424 


3.409 


3.410 



separation, the non-zero elements of the dielectric tensor, including local-field 
effects and the contributions of the non-local pseudopotential (cf. Appendix 
C), are e xx =5.1, e yy =5.5, and e zz =2.2 at the smallest imaginary frequency 
cj=0.036 Hartree, in agreement with effective-medium theory. Without the 
contributions of the non-local pseudopotential the values are e xx =5.9, e yy =6.5, 
and e zz =2.8, which underlines their importance for our test system. Varying 
the thickness or the separation produces changes in the quasiparticle energies 
that are of similar magnitude as the errors from an inadequate treatment of 
the anisotropy. In order to be able to investigate the surface or isolated-film 
limit, it hence proved to be essential to take the anisotropy modifications into 
account [33]. 

4-1 Convergence with respect to Z max 

In Table 1 we report the convergence of the quasiparticle energies with respect 
to the maximum value of / used for the evaluation of W lr in real space according 
to Equation (24). It can be seen that already with Z max = 2 the results lie within 
our level of accuracy (~ 0.05 eV), and with Z max = 4 absolute convergence is 
reached. These results also indicate that previous approaches [22,23], which 
have treated the anisotropy at the level of / max = 2, have incorporated the 
most important aspects of the anisotropy since the terms from I > 2 yield 
only minor corrections. Nevertheless, as the computational cost for evaluating 
higher terms in Equation (24) is negligible, we use Z max = 6 in practice. 

4-2 Convergence behaviour with respect to k-points 

In Figure 1 we show the convergence with respect to the number of k-points in 
the direction perpendicular to the surface for the quasiparticle energy of the 
lowest conduction state. Other states exhibit a similar behaviour. It is obvious 
that the original isotropic averaging for the screened interaction, notably in 
the long-range part, leads to an unphysical linear increase in the quasiparticle 
energy. In contrast, the anisotropic treatment converges rapidly. 
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Fig. 1. Convergence of the lowest conduction-band energy with respect to the num- 
ber of k-points N z perpendicular to the surface for (a) the original isotropic im- 
plementation and (b) with the anisotropy taken into account. The quantitative 
behaviour depends on the sampling in the parallel direction (N x = N y ). Note the 
different scales of the two graphs. 

The reason for the linear increase in the isotropic treatment is found in the 
inadequate treatment of the singularity, which is not fully removed. Integrating 
l/|k| 2 numerically yields for k x = k y = with Ak z = k max /N z 



n=l 



1 



{nAk z f 



dk 



1 



1 



1 



Ak z 



k 2 AK 



k„ 



(38) 



and hence a linearly diverging contribution, whose weight is proportional to 
Ak x Ak y ~ (N x Ny)~ l . When the k-sampling is increased in all three direc- 
tions simultaneously, no such linear divergence occurs, as can be seen from 
the three-dimensional plot in Figure 2 when going from the left side (small 
N x , N y ,N z ) to the right. However, such a restriction is undesirable and ineffl- 
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Table 2 

Extrapolated Ak x = Ak y — > quasiparticle energies in eV for the lowest conduction 
state for the isotropic and anisotropic treatment with a fixed number of k-points in 
the z-direction N z . 



N z 


1 


2 


3 


4 


5 


6 


isotropic 


3.063 


3.094 


3.090 


3.082 


3.074 


3.064 


anisotropic 


3.179 


3.160 


3.158 


3.157 


3.157 


3.156 



cient in practice. Therefore, only the proper anisotropic treatment enables us 
to investigate the importance of the k-point sampling in the direction perpen- 
dicular to the surface, which is directly related to the interaction with adja- 
cent slabs in GW calculations [33]. To our knowledge, the convergence in the 
perpendicular direction has not been addressed in previous GW calculations 
for slab systems, probably under the erroneous assumption that neighbouring 
slabs do not interact. 



The computed quasiparticle energies appear to be a linear function of the 
product Ak x Ak y when N z is kept fixed (not shown). This can be exploited 
to extrapolate the value for Ak x = Ak y — > 0, as shown in table 2. The ex- 
trapolated values evidently depend much less on the number of k-points in 
the ^-direction than those for finite Ak x and Ak y . For isotropic averaging, the 
situation greatly improves after extrapolation, but a small, systematic trend 
towards lower energies for larger N z is still present. This indicates that even 
after extrapolation no reliable convergence with respect to the number of k- 
points can be reached in the isotropic case. A comparison with the anisotropic 
treatment shows that the absolute error of the isotropic averaging is about 
0.1 eV, larger than could have been estimated from the isotropic data alone. 
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5 Summary 



We have presented a comprehensive account of the treatment of anisotropic 
screening in GW calculations that employ reciprocal space for the computa- 
tion of the screened interaction. In particular, we have demonstrated that this 
requires only small modifications of the original GW space-time implemen- 
tation [20]. The additional terms are computationally not very demanding. 
Furthermore, we have shown that the treatment of the anisotropy in other 
GW implementations can be understood in terms of approximations to the 
exact equations. 

The improvements presented in this article greatly increase the efficiency of 
GW calculations for anisotropic systems in the space-time method, e.g., for 
films and surfaces. This is mostly due to the fact that the fully anisotropic 
treatment enables us to converge the k-point sampling in the perpendicular 
and parallel direction separately, whereas isotropic averaging leads to an un- 
acceptable linear divergence in this case. The number of k-points required for 
converged results is thus reduced considerably. 
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A Angular expansion of vector expressions 

In this section we briefly summarise how simple expressions for a normalised 
vector k can be written in terms of spherical harmonics of the corresponding 
spatial angle Q^. The expansion of a scalar product k • r requires spherical 
harmonics of order / = 1 



i 




(A.l) 



m 



.=-1 
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with 




' o = \/ ^r z , r±i = Jy(F, + ir„) . (A.2) 



It is also straightforward to show that a tensor expression k T Lk, where L is 
symmetric, can be written in terms of spherical harmonics up to 1=2 as 

i 

k T Lk = E E ^m(fik) (A.3) 

le{0,2} m=-l 



with the coefficients 

-^00 — \ ~Tr{L X x + ^jm + L zz ) , L20 = \/ -r^{2L zz — L xx — , 

V 9 V 45 (A 4) 

/ 37T / 2"7T 

V to V to 



In the GW space-time method the term k T Lk appears in the denominator 
of the head element of the inverse dielectric matrix (9) as well as its product 
with an angular weight function iu(fi k ). In order to expand these expressions 
in spherical harmonics 



1 00 1 

rwr = E E ^lm^m(nk), (A.5) 

w(a) -£ E ^im(n k ), (a.6) 



k T Lk 



Z=0 m=-Z 



we determine the coefficients numerically by performing the following integrals 
on a Lebedev-Laikov grid [34] 



H lm = | d Q k ^(Q k )^-, (A.7) 

^, m = /dn k i^(n k )|^. (A.8) 

Since k T Lk and u>(O k ) are even functions, only even /-components contribute 
to the sums. 
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B Dielectric matrix 



For calculating the long-range limit of the symmetrised dielectric matrix we 
follow the derivation of Baroni and Resta [35]. For k — > we have 

Poo (k) ~|A; | 2 "head", (B.l) 
P G0 (k)~|A;| "wing", (B.2) 

P GG '(k) "body" . (B.3) 

This behaviour holds even for the exact polarisability of a non-metallic system 
and cancels the Coulomb singularity in the symmetrised dielectric matrix. In 
the context of the GW approximation Equations (B.l) to (B.3) are also valid 
for full self-consistency. However, here we restrict ourselves to the non-self- 
consistent case and derive expressions for the corresponding Taylor coefficients 
from the Adler-Wiser formula [36,37] for the polarisability 



W) - g/ A J^-S; (b.4) 



e i(k+G')-r 



where the sum over v and c runs over occupied and unoccupied states, respec- 
tively. For G' = and k — > this leads to 



«->-^S/ di '(^ (B ' 5) 

x (y?i>q|e~ lG ' r |y? cq ) (k • (<Poq|r|<P«q)) ' 
while for G = G' = and k — > the result is 

x (k • ((p m \r\(p CCL )) (k • ((p CCL \r\(p vq )) . 

The computation of the matrix elements (*|r|*) for the Kohn-Sham eigen- 
functions is presented in Appendix C. Writing the scalar products k ■ (*|r|*) 
as J2a k a (*| r a|*)> where a runs over the spatial (Cartesian) directions, we ar- 
rive at the following expressions for the wings of the symmetrised dielectric 
matrix in terms of a new vector quantity U(G,u;): 

^-gjb^s/ ^x^ (R7) 
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X 

k. 



-iG r 



V^cq) (VcqkalVuq) 



(B.8) 



and analogously for the head element in terms of a new tensor quantity F(u) 



e 00 (k, iu) 



^9 l k l (Svr)"^ 

rq 



^cq ^dc 



^cq ^vc 



a,/3 



k a k/3 



(B.9) 



(B.10) 



Since for many systems head and wings converge much slower with respect 
to the k-point sampling, but faster with respect to the number of conduction 
bands compared to the body [35], we compute them in a separate calculation. 
Equations (B.8) and (B.10) hold also for the self-consistent case, but the co- 
efficients U a (G;uj) and F a p{u) would have to be computed differently. The 
considerations following from now on then apply to self-consistent GW, too. 

Equations (B.8) and (B.10) illustrate that the limit for k — > is finite but 
will, in general, depend on the direction k/|k| = k in which the T-point is 
approached. We denote this directional dependence in the limit k = by 
the directional (spatial) angle f2 k - In existing implementations the treatment 
of the directional dependence varies. Sometimes, the direction f2 k is simply 
fixed to a single value. It is also common to carry the directional dependence 
through the inversion by performing a block- wise inversion [28], which is also 
the approach taken in the original space-time implementation [20]. For brevity, 
we omit the frequency argument from the following derivation. We denote the 
body of the symmetrised dielectric matrix (G ^ 0, G' ^ 0) at k = by B, 
the wings by w (a column vector) and w' 1 " (a row vector), and the head by H. 
The symmetrised dielectric matrix hence takes the form 



H(Q k ) w t(fi k ) 



v w(Q k ) 



B 



(B.11) 



Head, wings, and body of the symmetrised inverse dielectric matrix are then 
given by [28] 



^oo 1 (^k)=[^(^k)- E w G {Q k )B GGI w G/ (Q k ) 

^Go(^k) = -^oo(^k) E B GG' W G'(^k) , 



(B.12) 

(B.13) 
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^GG'(^k) — -BgG' + e Oo(^k)[ £ B GG" W G"(^) 



><[ E W G"( fi k)fi G ''G' 
G'YO 



Using equation (B.8), we now define the auxiliary vector 
S a (G) = E B^ G ,U a (G') 



(B.14) 



(B.15) 



and rewrite equation (B.12) as 

a,/3 G^O 



" k T Lk ' 



(B.16) 



thus defining L. Correspondingly, we have 



^Go(^k) - 



k-S(G) 



k T Lk 



£ GG' (^k) — Sqq, + 



k-S(G) 



k-S*(G') 



k T Lk 



(B.17) 
(B.18) 



C Kleinman-Bylander correction to the matrix elements of the 
position operator 



The matrix elements of r, which enter the expressions in the previous sec- 
tion, are in practice calculated via the commutator of r with the Kohn-Sham 
Hamiltonian h KS as 

Wr|^> = <M^il^S> . (C.1) 



While the contribution from the kinetic-energy operator is trivial to compute, 
that from the non-local pseudopotential V n \ is more cumbersome and has often 
been neglected in earlier calculations. We show here that it is possible to 
compute it efficiently in a separable expression. 

In its separable Kleinman-Bylander form [26] the non-local pseudopotential 
operator is written in the Dirac notation as 

Ki = £|<M-^<0,|, (C.2) 
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where fj, is a collective index {R M , l^, that runs over all pseudopotential 
projectors while M is in general given in a radial basis around a certain atomic 
position R M , i.e., 

M*) = fn»lM* ~ R„|) W^-R M ) • (C.3) 



In addition \i can run over chemical species, which does not alter the following 
derivation, except that /„; then also depends on the species. We will now 
show that the commutator can be factorised, which reduces the scaling to be 
linear in the number of plane waves instead of quadratic as demonstrated in 
a previous approach [38] . To this end we consider the commutator of r with a 
single projector 



(l^>4-<<^|r) - (r|^>4-(^l) 
^-[|^M„|(r-R„) - (r-R^I^X^I 



(C.4) 



Next we make use of the fact that r — R M can be expressed in the same radial 
basis as M 



[ r — R/Ja — | r ~~ R/i| E C am^m(^r-Rj , 



(C.5) 



m=— 1 



where a G {x, y, z) are the spatial directions, and c am yield the spatial com- 
ponents of the spherical harmonics for 1 = 1: 



Cam 


a = x 


a = 


m = 


-1 


i 


i 


m = 











m = 


1 


1 

V2 


i 
V2 





1 



We can then write the product in the radial basis, too, 



|0P := [r - R,U<M 
l 

= l r - R /J E c am^m(^r-R f ,)/n # .^(|r-R A1 |)lj f , n , M (fip-R M ) 
m=— 1 

= E E cS,^/^(|r-R M |)y iM (fi r _ R J (C.6) 

L=|Z M -1| M=-L 
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with 




^2c am >(lm \m\LM) , 



(C.7) 
(C.8) 



where (Im lm'\LM) is a Clebsch-Gordan coefficient. It is convenient to express 
0° in a plane- wave basis similar to what is done for M . If the functions f n i 
are given on a radial grid [39], f r nl is trivial to compute, and the same routines 
that are used to compute M (k + G) in the DFT calculation can be employed 
for the summands in 0^(k+ G). It must be emphasised that the sums over L 
and M contain only a very small number of non-zero terms (at most six). 

The final formula is thus again a separable expression 



The computational effort to set up a full N v x N c matrix for all three directions 
requires ANqN h(N v + jV c ) operations to calculate the ((f>n\(p v / c ) and ((f>u\(p v /c) 
projections and GN^NyNc operations to build up the 3 matrices from the 
projections in Equation (C.9). The scaling is thus linear in the number of 
G- vectors N G and not quadratic [38]. 

We have tested the size of the contributions from the non-local part of the 
pseudopotential for GaN and the II- VI compounds ZnO, ZnS, and CdS. We 
found that the macroscopic dielectric constant changes between +8 and - 
15%, which results in changes of the quasiparticle energies between -0.03 eV 
and 0.15 eV [3]. 
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